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ABSTRACT 

The dynamics of accretion discs around galactic and extragalactic black holes may 
be influenced by their magnetic field. In this paper we generalise the fully relativistic 
theory of stationary axisymmctric tori in Kerr metric of Abramowicz et al. l)1978|) by 
including strong toroidal magnetic field and construct analytic solutions for barotropic 
tori with constant angular momentum. This development is particularly important for 
the general relativistic computational magnetohydrodynamics that suffers from the 
lack of exact analytic solutions that are needed to test computer codes. 
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1 INTRODUCTION 

Accretion discs of black holes have been a subject of inten- 
sive observational and theoretical studies for decades. In 
most of the studies the discs are treated as purely fluid flows 
( Shakura , 



IFishbone fc Moncrief "19761 



Sunyaev 1974| 

lAbramowicz et al. 19781 IRees et al. 1982H . That is their 
magnetic fields are considered as dynamically weak and 
not important for the force balance determining the 
. disc structure. On the other, it has been long recog- 
nised that MHD-turbulence may hold the key to the 
nature of the disc viscosity that enables the inflow of 
matter into the black hole QShakura fe Sunyaev 1974| l. 
This expectation was confirmed with the discovery of the 
magneto-rotational instability jBalbus fe Hawley 1991| > 
that has become now one of the major topics in theoretical 
and numerical studies of accretion discs. In the simula- 
tions of radiatively inefficient discs they remain relatively 
weakly magnetized with the magnetization parameter 
f3 — (gas pressure) /(magnetic pressure) ~ 10 — 100 within 
the main body of the disc and rising up to j3 — 1 only near 
its inner edge and in the plunging region HHirose et al. 20041 
|McKinney fc Gammie 2004} . However, cooling accretion 
flows may well evolve towards the state with (3 < 1 where 
the vertical balance against the gravitational force of the 
central object is achieved by means of magnetic pressure 
alone IMachida et al. 20061 IPariev et al. 200311 . 

On the other hand, there is a general consensus 
that dynamically strong, ordered, poloidal magnetic field 
is required both for acceleration and collimation of the 
relativistic jets that are often produced by the astrophysical 



black hole-accretion disc systems. Such magnetic field 
is a key ingredient both in the models where the jet is 
powered by the black hole l |Blandford fc Znajek 19TT] > and 
in the models where it is powered by the accretion disc 
( Bisnovatyi-Kogan & Ruzma ikin| |Blandford fc Payne | l. 
The origin of this field is not very clear. The most popular 
idea is that it is carried into central parts of the disc-hole 
system by the accreting flow itself and that the net mag- 
netic flux gradually builds up there during the long-term 
evolution of the system jBisnovatyi-Kogan fc"R uzmaikin 
IThorne fc' Macdon alTJ 1982i . This" is more or less what 
is observed in recent numerical simulations of ra- 
diatively inefficient discs with initially weak poloidal 
field, (IHirose et al. 20041 |McKinney fc Gammie "2004| l . 
However, the strength of magnetic field that can be 
reached in this way is still unclear, e.g. ILivio et al. 19991 
[Uzdensky fc Spruit 2005] IMeier 20051 |McKinney 20051 . 
The most extreme case, where the disc pressure is dom- 
inated by ordered poloidal magnetic field is argued by 
Meier f2005$ as a model for the low/hard state of X-ray 
binaries. Although much weaker magnetic field is required 
to explain the observed power of quasar jets within the 
Blandford-Znajek theory, the pressure of poloidal magnetic 
field still has to be of the same order as the radiative pres- 
sure at the inner edge of the radiatively supported accretion 
discs (Bcgclrna n et al. 1984| . Bogovalov and Kel'ner (20051 
constructed a model of accretion disc where its angular 
momentum is carried away by a magnetized wind alone and 
the inflow of matter and magnetic field is driven entirely by 
the magnetic torque applied to the disc by the wind. These 
and other results show that the structure of accretion discs 
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with dynamically strong ordered magnetic is a matter of 
significant astrophysical interest. 

General steady-state axisymmetric magnetized dicks 
are described by a very complex equation of the Grad- 
Shafranov type which does not allow general analytic 
solutions ILovelace et al. 19861 . However, the equilibrium 
models with pure azimuthal (toroidal) field are much 
simpler and can be constructed using the approach de- 
veloped for unmagnetised discs (Abramowicz et al. 19781 
IKozlowski et al. 197811 . In fact, the predominant motion in 
accretion discs is a differential rotation and one would ex- 
pect the azimuthal magnetic field to dominate in the interior 
of such discs. 

Okada et al.(1989l constructed a particular exact so- 
lution for the problem of equilibrium torus with purely az- 
imuthal magnetic field. In their approach the Paczynski & 
Wiita 119801 potential is utilised to introduce the black 
hole gravity and the flow dynamics is described by the non- 
relativistic equations. Although this approach seems to work 
fine in the case of non-rotating black holes, it is no longer 
satisfactory in the astrophysically important case of rapidly 
rotating black hole where full relativistic treatment is re- 
quired. In this paper we develop the general relativistic the- 
ory of magnetized tori around rotating black holes. Our anal- 
ysis closely follows the one of Kozlowski et al. 119781 and 
Abramowicz et al. 119781 . including the notation. 

Even if tori with dynamically strong magnetic fields 
do not exist in nature the exact solutions presented here 
can still be useful for numerical general relativistic magne- 
tohydrodynamics, GRMHD, which has attracted a lot of 
interest recently fKoide et al. 19991 Komissaro v 200ll 

De Villiers fc Hawley 2003] IGammie et al. 20031 

Komissarov 20041 IDuez et al. 20051 lAnton et al. 200611 . 
One of the problems with general relativistic computa- 
tional magnetohydrodynamics is the lack of exact analytic 
and semi-analytic solutions that could be used to verify 
computer codes. The current list of such solutions includes 
1) the spherical accretion onto a nonrotating black hole 
with monopole magnetic field where the magnetic field is 
dynamically passive IKoide et al.1 999 1 , 2) the perturbative 
force-free solution of Blandford-Znajek 119771 . where both 
the pressure and the inertia of matter is neglected and the 
black hole rotates slowly, and 3) the so-called Gammie flow, 
that is a one- dimensional solution of the Weber-Davis type 
that applies only to the equatorial plane of a black hole 
IGammie 19 99). One may also try to utilize the axisymmet- 
ric semi-analytic self-similar solutions constructed recently 
by Meliani et al. 120051 but those apply only along the 
symmetry axis. Obviously, the magnetized torus solution is 
very welcome an addition to this scarce list. In fact, it is 
the only solution that is not only a) fully multi-dimensional 
and b) involves dynamically important magnetic field, but 
c) also applies in the case of a rapidly rotating black hole. 
In fact, we have already used this solution to test the 2D 
GRMHD codes described in Komissarov (EUtJTl ISTMll . 

Throughout the paper we use the relativistic units 

where c = G = M = 1 and ( h ++) signature for the 

space-time geometry. Following Anile (119891 the magnetic 
field is rescaled in such a way that the factor 4ir disappears 
from the equations of relativistic MHD. 



2 BASIC EQUATIONS 

The covariant equations of ideal relativistic MHD are 

V Q T Q ' 3 = 0, 
V Q *F aS = 0, 
V a p« a = 
where 



V = (w + b 2 )u a u + (p + ^b 2 



a/3 



<xiB 



b a b 



(1) 

(2) 
(3) 

(4) 



is the energy- momentum tensor, w, p and u a are the fluid 
enthalpy, pressure and 4-velocity of plasma respectively, and 
g a /3 is the metric tensor IDixon 1 978 Anile 1989J . 

* F af3 = h c u _ b fS ua (g-) 

is the Faraday tensor and b a is the 4-vector of magnetic 
field. In the fluid frame b a = (0,B), where B is the the 
usual 3- vector of magnetic field as measured in this frame, 
and thus 



u a b a — 0. 

In the following we assume that 



(6) 



(i) the space-time is described by the Kerr metric and 
that {t, 4>, r, 9} are either Kerr-Schild or Boyer-Lindquist co- 
ordinates, so that 

SW = g^,4> = 0; (7) 

(ii) the flow is both stationary and axisymmetric, so that 

/,* = U = (8) 

for any physical parameter /; 

(iii) the flow is a pure rotation around the black hole, that 

is 

u r = u B = 0; (9) 

(iv) the magnetic field is purely azimuthal, that is 



0. 



(10) 



In terms of partial derivatives the continuity equation 
reads 

where g is the determinant of the metric tensor. Given the 
symmetry conditions J7J and @ this equation reduces to 



(n 



-gpu ),j 



0, 



where here and throughout the whole paper i = r, 9. Fi- 
nally, the condition tells us that this equation is always 
satisfied. 

Since the Faraday tensor is antisymmetric one may re- 
duce the Faraday equation to 

It is easy to see that this equation is also automatically 
satisfied given the conditions 171101 . 

Thus, the only non-trivial results follow from the 
energy-momentum equation Contracting this equation 
with the projection tensor h a g = S a p + u a u/3 we obtain 



(w + b 2 )u v u ^ + (p + & 2 ),i — Kb"^ = 0. 



(11) 
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In terms of the angular velocity, 



and the specific angular momentum 

I = -u^jut. 



(12) 



(13) 



(both ut and u$ are constants of geodesic motion) this equa- 
tion reads 



(In \u t \ 



n 



HI w ILw 



(14) 



where 

£(r, 6) = g t 4,gt<i> - gttg<M>- (15) 

When b 2 — > equation 114H reduces to equation 7 in 
Abramowicz et al. 11978(1 which describes equilibrium non- 
magnetic tori. Thus, the Lorentz force vanishes and we have 
a force-free magnetic torus provided 

Lb 2 = const. 

Far away from the black hole <?t</> — * 0, g^ — > r 2 sin 2 8, 
gtt—*—l an d this equation reduces to the familiar 

4 const 



r sin 6 



p = Pm. = 0, Ut=u tin , l = lin 

one finds the constant of integration as 

hn 

const = In lilt . I — 



Qdl 



I -in' 



(22) 



(23) 



Following, Abramowicz et al. ({19780 we introduce the total 
potential, W, via 



W = In I ii t I 



(24) 



where loo is the angular momentum at infinity. Provided 
that loo is finite we have u too = — 1 and Woo = 0. Using the 
total potential we can rewrite eg. 12111 as 



w-w in + I d -P+ [<^=0. 



W J It) 

o o 



(25) 



With exception for the last term this equation is the same 
as eq.9 in Abramowicz et al.(1978) and eq.24 in Kozlowski 
et al.iT575t. 



3 INTEGRABILITY CONDITIONS 

For a barotropic equation of state, where 
w — w(p), 

equation l|14|l leads to 



Q, M _ d{Cb 2 ) 



^ln \ut\- 










o 1 



2Cw 



(16) 



(17) 



In the case of a non-magnetic torus this equations implies 
that 



Q, = fl(l) 



(18) 



and thus the surfaces of equal 57, I, p, and p co- 
incide lAbramowic z 197T1 lAlpramowic z et al. 1978^ Obvi- 
ously, this does not have to be the case for magnetized tori. 
If, however, we still assume that fi = fi(Z) then ea. H7H can 
be written as 



d I In lit | + 



w J 1 





Qdl 



d(Cb 2 ) 
2Cw 



(19) 



that implies that the expression of the right-hand side is a 
total differential. Hence, 



W = U>(p~m), 



(20) 



where iD = Cw and p m = Cp m , where p m = b 2 /2 is the 
magnetic pressure, and ea. l|19^ integrates to give 



lu|«,|- / 



ndi 

w j i-m 





dp m 
w 



= const 



(21) 



Assuming that on the surface of the disc, and hence on 
its inner edge, 



4 BAROTROPIC TORI WITH CONSTANT 
ANGULAR MOMENTUM 

4.1 Theory 

Here we adopt particular relationships w = w(p), Q — O(Z), 
and w — w(pm) that allow to express the integrals of equa- 
tion 1251 in terms of elementary functions. Namely, we as- 
sume that 



I = la, 
p = Kw K 

Pm = K m W V . 

The last equation can also be written as 

p m = K m C V ~ 1 W V . 

Then eq. l25H reduces to 



W - Win + 

where 

W = \n\utl 



ft — 1 it) 



J! Pm 

— 1 It) 



0, 



(26) 
(27) 
(28) 

(29) 
(30) 
(31) 



The obvious parameters of the model are k, r/, lo, and 
Win- Two more parameters are needed. We chose these to 
be the enthalpy, w c , and the magnetization parameter f3 c — 
(p/Pm)c at the disc centre, r = r c . Following Abramowicz et 
al.( 1978 1 we define the disc centre as one of the two points in 
the equatorial plane where lo equals to the Keplerian angular 
momentum 



±(r 2 =F 2ar 1/2 + a 2 
r 3/2 _ 2r i/2 ± a 



(32) 



where the upper sign is used if lo > and the lower sign 
otherwise (Ba rdeen et al. 197211 . The second point is to the 
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disc "cusp" , r cusp < r c ( Abramowicz ct al. 1978 1 . There ex- 
ist a number of obvious constraints on the values of these 
parameters. 

In order to avoid divergence of the second and/or the 
fourth terms in 130H at the disc surface one should have 

k,/3 > 1. 

Under this condition the disc surface is fully determined by 
the choice of Win and does not depend on the disc magneti- 
zation. This property has already been noticed in Okada et 
al. 115591 . 
Only if 

|Zo| > \lms\, 

where l ms is the radius of the marginally stable Keplerian 
orbit CBardee n et al. 1972L the disc is detached from the 
event horizon (Abramowicz et al. ( 1978). Solutions attached 
to the event horizon are improper because they diverge at 
the event horizon and cannot be continued through it.) The 
value of lo determines the total potential which can be writ- 
ten as 



W(r,6) = - In 
where 



and 



(33) 



A = g^ + 2logt<p + ?o5« 

(see ea lXTol l. 

If |2o| > \lmb\ then the disc has finite outer radius only 

for 

Win < 0. 

If \lms\ < | Jo | < \lmb\ where l mb is the angu- 
lar momentum of the marginally bound Keplerian orbit 
(Bardccn ct al. 19721, then the disc remains detached from 
the black hole only if 

Win < Wcusp, 

where W cusp is the value of the total potential at the cusp 
ramowicz et al. 197811 . 
From eg. (1301 one finds the gas pressure 



Pc = w c (W ln - W c ) ( — ^— + 

r/-l /3 c 

and then the magnetic pressure 

Pm c = Pc/Pc 



(34) 



(35) 



at the disc centre. Using these one finds the constants K 
and K m of the barotropics I|27I29[I . 

Now one can compute the solution at any location in- 
side the disc. Given the coordinates (r, 6) one computes the 
potential W using equation l|33|l . If W < Win then this 
point is inside the disc. The next step is to find the en- 
thalpy, id, as a solution of eg. 1301 . Once w is found one can 
find p and p m from eos. l|27^ and 129|l . The 4- velocity vector, 
u v = (it*, u^, 0, 0), is given by 

— fht 4 , 



Model 


k 


Tcusp 


r c 


W^cusp 


W c 


w m 


/3c 


A 


2.8 


1.58 


4.62 


0.702 


-0.103 


-0.030 


0.1 


B 


2.6 


1.78 


3.40 


-0.053 


-0.136 


-0.053 


1.0 



Table 1. Parameters of the models used for test simulations 



where Q, can be found via 



gt4> + gulp 
g<t>4> + gt<i>lo 



(see Appendix.) The 4-vector of magnetic field, b v 
(&*, b*, 0,0), is given by 



b* = ±y / 2p m /A, b l = lob* 



(see eq |A16t . 



u t (l — loty 



4.2 Simulations 

In order to illustrate the usefulness of this solution for test- 
ing GRMHD computer codes we constructed two equilib- 
rium models and used them to setup the initial solution for 
2D axisymmetric simulations using the code described in 
Komissarov 121) (J4t . In both models the black hole has spe- 
cific angular momentum a = 0.9 which gives r m b = 1.73, 
r ma = 2.32, l m b — 2.63, l ms = 2.49. In both cases, we used 
the same value for the barotropic powers «, rj — 4/3 and the 
polytropic equation of state with the same ratio of specific 
heats, 7 = 4/3. The other parameters are given in Table 1. 
Notice, that lo > l m b in model A, whereas in model B one 
has lms < lo < lmb and this allows accretion through the 
torus cusp ( Abramowicz et al. 19781 IKozlowski et al. 197811 . 

Initially, the space outside of the tori is filled with a 
rarefied non-magnetic plasma accreting into the black hole. 
Its density and pressure are 

p = 10 _3 /9 C exp(— 3r/r c ), p = Kp K . 

Its velocity in the frame of local fiducial observer, FIDO, 
frh orne fc Macdonald 1982 ) is radial and has the magni- 
tude 

« = ^(l-(r s /r) 4 ) 

where (3 r is the radial component of velocity of the spacial 
grid relative to FIDO. This introduces inflow through the 
horizon with the local velocity of FIDOs and at the same 
time allows only small poloidal velocity jump at the torus 
surface. The latter property allows to avoid strong rarefac- 
tions that may originate at the torus surface right at the 
start of simulations. The initial distributions of p, £1, and /3 
for the model A, that is a more extreme case, are shown in 
figure Model B appears very similar. 

In these simulations we utilized the Kerr-Schild coordi- 
nates (r, 6). The computational domain is [1.35, 53.3] x [0, 7r] 
with 320 cells in each direction. At lower resolution the ef- 
fects of numerical diffusion become rather noticeable. The 
grid is uniform in the ^-direction and the cell size in the ra- 
dial direction is such that in the equatorial plane geeA8 2 — 
g rr Ar 2 . This ensures that computational cells have approx- 
imately equal lengths in both directions and throughout the 
whole grid. 

The rotational period of the disc centre is rather long 
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Figure 1. Exact solution for model A. The horizontal coordinate is R = rsinff and the vertical coordinate is Z = rcosS Left panel: f3, 
Middle panel: Q, Right panel: p. 




Figure 2. The density distributions for the model A obtained with our actual numerical code and its corrupted versions. In all these plots 
the linear scaling is used. Left panel: the solution obtained with the uncorrupted code at t = 200. Middle panel: the solution obtained 
with the corrupted version of the code at t = 14. Here we retained only the contribution of magnetic pressure in the Maxwell stress 
tensor. Right panel: the solution obtained with the corrupted version of the code at t = 14. Here all the electromagnetic contributions to 
the stress-energy-momentum tensor were omitted. 



r r — 68 in model A and r r = 45 in model B. However, 
the dynamical timescale, Td, for the disc centre is signif- 
icantly shorter. Indeed, in both of these models the fast 
magnetosonic speed in the disc centre is aj ~ 0.18, whereas 
the length scale of pressure distribution at this location is 
L a — 1.0. This gives the dynamical timescale Td — 6. Since it 
is the dynamical timescale that determines how quickly the 
system reacts to perturbations of its equilibrium state, it is 
quite sufficient to carry out simulations for only t = few Td in 
order to test the ability of our numerical code to correctly 
reproduce these equilibrium solutions. In order to demon- 
strate this we have carried out test simulations for model A 
not only with our "proper code" but also with its two cor- 
rupted versions. In the first version, it is only the electromag- 
netic pressure, (B 2 + E 2 )/2, that was taken into account in 
the calculations of the Maxwell stress tensor. In the second 



and more drastic version we omitted all the contributions 
of the electromagnetic field to the stress-energy-momentum 
tensor. The results are presented in figure One can see 
that already at t ~ 2rd the solutions obtained with the cor- 
rupted versions are significantly different from the initial 
equilibrium solution. With the second corrupted version the 
disc simply collapses towards the equator due to the lack of 
magnetic support against gravity. 

Figure [5] also presents the proper numerical solution at 
t=200, which is 30 -f- 40 times larger than Td- "Naked eye" 
inspection shows that this solution is very similar to the 
initial equilibrium one that is presented in figure This is 
confirmed by the ID density plots showing in details the 
distributions along and across the symmetry plane: figure |3] 
for model A and figure [I] for model B. Other images also 
reveal what appears to be surface waves propagating away 
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Figure 3. Density distribution at t = 60 for model A. The solid lines show the exact equilibrium solution. Left panel: log 10 p against 
polar angle 9 at r = r c , the disc centre location. Right panel: p against r in the equatorial plane. 





Figure 4. Density distribution at t = 60 for model B. The solid lines show the exact equilibrium solution. Left panel: log 10 p against the 
polar angle 8 at r = r c , the disc centre location. Right panel: p against r in the equatorial plane. 



from the black hole and becoming noticeable in the more 
remote parts of the tori. 



5 SUMMARY 

In this paper we have generalized the relativistic 
theory of thick accretion discs around Kerr black 
holes (IFishbone fc Moncrief 19761 lAbramowicz et al. 19781 
Kozlowski et al. 197811 by including dynamically strong 
toroidal (azimuthal) magnetic field. As expected, this in- 
clusion of magnetic field leads to a whole new class of equi- 
librium solutions that differ in strength and spacial distri- 
bution of the field. In particular, we have described the way 
of constructing barotropic tori with constant angular mo- 
mentum - under such conditions the differential equations 
of magnetostatics are easily integrated and reduce to simple 
algebraic equations. 



By now it has become clear that magnetic field plays 
many important roles in the dynamics of astrophysical black 
hole-accretion disc systems. However, the structure of this 
field may be more complex than just azimuthal loops and 
as the result the analytic solutions constructed in this paper 
may turned out not to be particularly suitable for mod- 
elling real astrophysical objects. Even so they will be cer- 
tainly helpful in testing computer codes for general rela- 
tivistic MHD that are becoming invaluable tools in the as- 
trophysics of black holes/accretion discs. 
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APPENDIX A: 

Here we derive equations 1111 and 1 1 H of the main paper. 
Al Equation lfTT]l 

The contraction of the energy-momentum equation, 

V Q T Q/3 = 0, 
with the projection tensor 

t a ca i a 

f3 — 13 + u U P 

leads to 
/i Q t V 7 Uw + b 2 )u' Ua \ + 

^V 7 ((p+6 2 /2)5 7 Q ) + (Al) 
h a l V 1 (6 7 6 a ) = 0, 

where i — r,9. Applying the product rule and using that 
h a fjU a — and u a V 7 u a = one can write 



/i"V 7 ((iu + b 2 )u 1 u a ) = (w + b 2 ) [V 7 (w 7 Wj) — Ui\7ju"'] . 

The second term in the square brackets vanishes because of 
the symmetry conditions 171101 . Applying the well known 
result for the divergence of a second rank symmetric tensor 
on obtains 



v^g 2 

The first term in this equation vanishes because of the same 
symmetry conditions whereas 

it V n V 

g^,iwu = — iu v u ti . 

Thus, one has 

/i°iV 7 ((w + 6 2 )u 7 u a J = (w + b 2 )u v u / ^ i . (A2) 
Now we deal with the second term in 1AH . 
ft Q 4 V 7 ((p + b 2 /2)P a ) = / I 7 V 7 (p + 6 2 /2) 

= {p + b 2 /2), l +u l u<{p + b 2 /2), 1 . 



The second term on the right-hand side of this equation 
vanishes because of the symmetries and hence we have 

|2 . 



h a ^((p + b 2 /2)P a ) =(p + b 2 /2), 
As to the third term in lAH . 

1 



(A3) 



/i Q iV 7 (lPb a ) = h° 



y & Q ), 7 — -g^,ab^b l 



{g^,ib^ l b v + UiU a # MI , lQ 6 M 6") 



-i(6 2 ,i-2W, i ) 



(A4) 



During this reduction we twice used the symmetry condi- 
tions. Substituting of results (IA2IA4I into eq. lAH leads to 



(w + b 2 )u v u \i + (p + 6 2 ),i — b„b" ti = 0, 
which is equation 1111) of the main paper. 



(A5) 



A2 Equation 111411 

From the definitions of the angular velocity, fi, and the an- 
gular momentum, I, and the symmetries of the problem it 
immediately follows that 

9t$ + g<t><j>Q 



I = 



gtt + fftr>^ ' 

gt<t, + gul 



(A6) 



(A7) 



94>4> + gt<i>l ' 

where we assumed the Kerr metric in either Kerr-Schild 
or Boyer-Lindquist coordinates. From the definition of 4- 
velocity it follows that 

II V 1 

g^vUu = — 1. 



In the case under consideration this leads to 

{uf = -{g tt + 2g ttj ,n + g^ti 2 )' 1 , 
and 



u u t 



l-in 

From HA7IA9|I one finds 
^2 C 



(my 



A' 



where 

£ — gt4>gt<t> — <7tt<?00i 

and 

A — g<j,<t> + 2lg t ^ + I 2 gtt- 



(A8) 



(A9) 



(A10) 



(All) 



(A12) 



From the constraint equation © and the conditions 19I1M 
one finds 



6* = lb*, 
b t = — Ofo^, 
These allow to write 
b 2 = b^il- ISl), 



and 



(6*) 2 A 



(A13) 
(A14) 

(A15) 

(A16) 



Now we reduce the terms u v u v ti and b^b 1 ' i in equation 1A5I . 
Using the definitions of £1 and I one derives 



= — u Ut,i — U W0,i 

= — U (Ut,i — n(lui),i) 

= — u* («t,i(i — in) — nuti,i) ■ 

The substitution of w* from IA9t into the last equation gives 



U„u" t i = (In \u t \),i — 



Q 



-I 



(A17) 



i -in 

Using t£T5l and llATHl one writes 
Kb" j = btb* :i + b^F^ 

= -nb^,(ib*),i +b^ i 

= 6 6^(l - in) - Slb^b%. 
The substitution of &</, from IA15I into this equation gives 
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Substituting ljA17IA18> into ljA"5jl one finds 

+ b 2 i = 
Using llATol and dXTBI one obtains 

,2 f, "t \ , 2 (£b 2 ),, 

Thus, fKTfy can be written as 

which is equation 114H of the main paper. 



(A18) 



(A19) 



(A20) 
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